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Abstract 


An  efficient,  stable,  explicit,  first-order  cross-stream  differencing  scheme  having  low  truncation 
error  is  derived  and  applied  to  three  dimensional  integral  boundary  layer  equations.  The 
analysis  is  considered  in  detail  for  the  particular  case  of  the  momentum  integral  equations  and 
Head’s  entrainment  equation  with  power  law  profiles  and  Mager’s  cross  flow  assumption.  In 
comparison  with  upwinding  differencing,  the  new  scheme  has  better  stability  properties,  smaller 
truncation  error,  and  smaller  artificial  viscosity. 


Resume 

Un  schema  de  differentiation  transversale  de  premier  ordre  explicite, 
stable  et  efficace,  avec  une  faible  erreur  de  troncature,  est  etabli  et 
applique  aux  Equations  integrales  tridimensionnelles  de  la  couche  limite. 
L* analyse  est  consideree  en  detail  dans  le  cas  particulier  des  Equations 
integrales  de  quantity  de  mouvement  et  de  1' Equation  d'entrainement  de 
Head,  avec  des  gradients  exponentiels  et  l'hypothese  d*6coulement 
transversal  de  Hager.  Compare  a  la  differentiation  longitudinale,  le 
nouveau  schema  possede  une  plus  grande  stability ,  une  plus  petite  erreur 
de  troncature  et  une  viscosite  artificielle  plus  faible. 
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2.1  Results  of  integration  with  /(0,y)  =  1  -  cos(jry) 

2.2  Results  of  integration  with  /( 0,y)  =  (1  —  |1  —  y|)/2 


4.1  Components  of  G  as  functions  of  H 


^ I 


Notation 


A,B  :  Matrices  of  a  hyperbolic  system  of  first  order  partial  differential  equa¬ 
tions. 

C  :  Matrix  defined  in  equation  (4.17). 

D  :  Matrix  defined  in  equation  (3.11). 

D+  ,D~  :  Forward  and  backward  differencing  operators  defined  by  equation  (2.4). 
E  :  Matrix  defined  in  equation  (4.17). 

Et  :  Truncation  error. 

/  :  Independent  variable  in  the  convection  equation. 

fk  :  Independent  variables  in  the  multi-dimensional  hyperbolic  system  of 
equations. 

:  Fourier  component  of  /  at  the  point  (xt,j /,■). 

F  :  Matrix  defined  in  equation  (4.16). 

G  :  Matrix  equal  to  EC-1. 

hii,hi2,hii,h22,ki,k2  :  Boundary  layer  velocity  profile  functions. 

H  :  Shape  factor. 

P  :  Stability  factor  definined  by  equation  (2.7). 

Ql,Q2,QLl  '•  Independent  variables  in  the  boundary  layer  equations. 

Qx2 ,  E1 ,  E2  :  Variables  of  which  cross-stream  derivatives  are  taken  in  the  boundary 
layer  equations. 

t  :  Tan  of  the  cross  flow  angle. 

T  :  Matrix  defined  in  equation  (4.16). 

v(n)  :  The  nth  eigenvector  of  A-1B. 

V *  :  Velocity  in  the  x  direction. 

Vv  :  Velocity  in  the  y  direction. 

Vl,V2  :  Contravariant  components  of  the  potential  flow  velocity. 
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V -L1,Vi2  :  Contravariant  components  of  a  vector  perpendicular  to  the  potential 
flow  velocity  but  having  the  same  magnitude. 

wM  :  Vectors  defined  by  equation  (3.5). 

x,  s  :  Non-orthogonal  surface  coordinates  in  which  the  boundary  layer  equa¬ 
tions  are  solved. 

x,y  :  Coordinates  in  which  the  convection  equation  is  solved. 

As,  Ax,  Ay  :  Grid  spacing  in  the  8,  x,  or  y  directions. 

a  :  Coefficient  defining  the  cross-stream  differencing  scheme. 

0  k  :  Coefficients  giving  /„  as  a  linear  combination  of 
6  :  Boundary  layer  thickness. 

6mn  :  Kronecker  delta. 

A„  :  The  nth  eigenvalue  of  A-1B. 

/i  :  A  parameter  defined  in  equation  (2.12). 

Bold  face  characters  are  reserved  for  use  as  vectors  or  matrices. 
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1  Introduction 


Numerical  integration  of  a  system  of  boundary  layer  equations  is  widely  used  for  the 
prediction  of  the  flow  over  ship  hulls  and  aircraft.  When  the  boundary  layer  equations  are 
written  in  integral  form  and  solved  in  a  coordinate  system  roughly  aligned  with  the  flow,  a 
popular  method  of  discretization  is  a  first  order  explicit  scheme  using  upwind  differencing  for 
the  cross-stream  derivatives.  This  method  has  the  virtue  of  extreme  simplicity;  however,  a 
disadvantage  of  the  method  is  that  it  is  unstable  when  the  direction  of  forward  integration 
lies  between  the  potential  flow  streamlines  and  the  limiting  streamlines  at  the  hull  surface. 
While  the  occurrence  of  instability  will  normally  be  infrequent  and  will  usually  not  persist,  it 
is  most  likely  to  occur  when  the  cross-flow  is  high:  that  is,  at  a  time  when  the  modelling  of 
the  boundary  layer  is  most  suspect.  When  testing  the  validity  of  a  boundary  layer  prediction 
method,  it  can  be  difficult  to  ascertain  whether  problems  occurring  at  regions  of  high  cross-flow 
are  due  to  deficiencies  in  the  boundary  layer  modelling  or  due  to  instabilities  in  the  integration 
procedure.  A  completely  stable  differencing  scheme  would  remove  this  ambiguity. 

A  second  disadvantage  of  upwind  differencing  is  that  it  is  discontinuous  with  respect  to 
the  direction  of  the  flow:  i.e.  an  infinitesimal  change  in  the  flow  which  causes  the  cross-stream 
component  of  the  flow  to  change  sign,  will  cause  finite  changes  in  the  solution  due  to  the  change 
in  the  direction  of  cross-stream  differencing.  If  the  boundary  layer  equations  are  being  used  in 
an  iterative  loop  to  account  for  the  interaction  between  the  boundary  layer  and  the  potential 
flow  (see,  for  example,  Chapter  8.3  in  Reference  [1]),  this  discontinuity  can  inhibit  convergence. 

In  this  memorandum  an  improved  scheme  for  the  cross-stream  derivatives  is  proposed. 
The  new  scheme  is  also  of  first  order,  is  stable  for  sufficiently  small  step  size,  and  is  continuous 
with  respect  to  the  direction  of  flow.  It  also  has  smaller  truncation  error  and  artificial  viscosity 
than  the  upwind  differencing  scheme. 

Motivation  for  the  new  differencing  scheme  for  the  boundary  layer  equations  is  given  by 
first  examining  the  pure  convection  equation  with  one  dependent  variable.  For  this  equation,  the 
new  scheme  is  shown  to  be  the  optimal  explicit  differencing  scheme  in  that  it  is  (conditionally) 
stable  and  minimizes  the  truncation  error.  The  scheme  is  then  extended  to  systems  of  equations 
(i.e.  equations  with  more  than  one  dependent  variable).  Finally,  use  of  the  scheme  in  the  context 
of  the  boundary  layer  equations  is  discussed.  For  this  purpose  the  boundary  layer  equations 
are  assumed  to  be  the  two  momentum  integral  equations  and  Head’s  entrainment  equation; 
these  are  the  equations  used  by  the  HLLFLO  computer  programs  developed  at  DREA  for  the 
prediction  of  the  flow  around  ship  hulls[l|.  A  complete  discussion  of  these  equations  is  given 
by  Hally (2], 


2  Stability  of  the  Convection  Equation 


Consider,  first,  the  simple  one  dimensional  convection  equation, 

vtj-+vv^r  =  ° 

ox  ay 


(2.1) 


where  Vx  will  be  assumed  positive  with  Vx  >  |V*,|.  To  solve  the  equation  numerically,  a 
rectangular  grid  (x*,yy)  is  used.  For  simplicity,  the  cross-stream  step  size  will  be  assumed 
uniform:  i.e.  y;+i  -  yy  =  Ay  for  every  j.  The  derivatives  are  approximated  by 


where  Ax  =  x*+i  —  x*  and 


|£  - 

fk+i,i  _  jk,j 

(2.2) 

dx 

Ax 

f  - 

dy 

\  V/+12V/ 

(2.3) 

fk,i+ 1 

A 

-  /*'■>■  D.f  fk'i  -  f'J-1 

(2.4) 

The  upwind  differencing  scheme  has  a  =  sgn(Vv).  Central  differencing,  which  has  accuracy  of 
second  order  in  the  y-derivatives,  has  a  =  0. 

The  stability  of  a  differencing  scheme  is  determined  by  supposing  a  solution  of  the  form 


fk,}  _  pkjj&v 

Substitution  into  equation  (2.1)  and  using  equations  (2.2)  and  (2.3)  yields 

fk+l _  pk,j 


1  -  §1(1  -  a)(«‘4*  -  1)  +  (1  +  «)(1  -  «-“")] 


where 


P  = 


V*Ax 

V*Ay 


The  scheme  is  stable  if  |F*+lj’|  <  |Fkj|  which  occurs  when 

aP  <  1  and  P( a  —  P)  >  0 


(2.5) 

(2.6) 

(2.7) 

(2.8) 


(see,  for  example,  Peyret  and  Taylor[3]).  Hence,  there  is  stability  only  if  the  scheme  is  weighted 
towards  upwind  differencing:  i.e.  sgn(a)  =  sgn(P).  Complete  upwind  differencing  is  stable  if 
|P|  <  1;  this  is  called  the  Courant-Friedrichs-Lewy  (CFL)  condition.  Central  differencing  is 
unstable. 


(2.9) 
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The  truncation  error  of  the  scheme  of  equations  (2.2)  and  (2.3)  is 


Et  =  V 


2  dx2 


-  aV* 


AydV 

2  3y2 


+  0(  Ax2,  Ay2) 


The  term  in  Ay  has  a  dissipative  effect  if  aVw  >  0.  Since  this  term  is  physically  similar  to 
a  diffusion  term,  it  is  often  called  “artificial  viscosity”.  For  a  first  order  scheme  the  artificial 
viscosity  is  necessary  for  stability,  but  to  avoid  unwanted  diffusive  effects,  it  should  be  as  small 
as  possible.  If  Vx  and  Vv  are  constant,  equations  (2.1)  and  (2.7)  can  be  used  to  rewrite 
equation  (2.9)  as 

V-Ay(a-P)  31/  ,  „,A  ,  .  2. 


E,= 


dxdy 


+  0(  Ax2,  Ay2) 


(2.10) 


In  order  to  minimize  the  truncation  error,  equation  (2.10)  suggests  that  one  use  a  =  P, 
a  scheme  which  will  be  termed  partial  upwinding  for  the  purposes  of  this  memorandum.  This 
scheme  is  stable  provided  that  |P|  <  1,  the  CFL  condition;  it  therefore  requires  no  smaller  step 
sizes  than  complete  upwind  differencing.  Moreover,  since  a  is  continuous  with  respect  to  P, 
the  differencing  scheme  is  continuous  with  respect  to  Vv .  The  artifical  viscosity  is  also  reduced 
in  comparison  with  that  of  the  complete  upwind  differencing  scheme. 

In  a  more  general  convection  equation,  A  and  B  may  not  be  constant  and  there  may  be 
production  terms. 

,  D  „«/ 


Akn(x,  y,  /)  —  +  Bkn(  X,  y,  /)  —  =  Ck{x,  y,  /) 


(2.11) 


The  differencing  scheme  defined  by  equations  (2.2),  (2.3),  and  a  =  P  is  stable  when  used  for 
equation  (2.11)  but  will  have  truncation  error  which  is  first  order  in  Ax  and  Ay.  However,  if 
convection  is  the  dominant  mechanism  describing  the  evolution  of  the  solution,  the  truncation 
error  will  be  small. 

The  relative  accuracy  of  central  differencing,  upwind  differencing,  and  partial  upwind 
differencing  have  been  compared  using 


V*_ 

V1 

The  solution  to  equation  (2.1)  is  then 


My(i  -  y) 

l  +  /ii(l-2y) 


/(*.  y)  =  /(°.y  +  ^xy{l  -  y)) 


(2.12) 


(2.13) 


In  this  flow,  Vv  is  zero  along  the  coordinate  lines  y  =  0  and  y  =  1;  hence,  specification  of  the 
value  of  /( 0,y)  for  0  <  y  <  1  is  sufficient  to  determine  the  value  of  /  for  any  x  >  0.  In  this 
sense,  the  flow  mimics  the  evolution  of  the  boundary  layer  on  a  ship  hull  which  is  confined 
between  streamlines  at  the  keel  and  at  the  waterline  (see  Hally  [2]). 

The  parameter  can  be  used  to  alter  the  relative  magnitude  of  the  components  of  the 
velocity  in  the  x  and  y  directions.  Since  one  cam  show  that 


A*y(i  -  y) 


1  +  nx{  1  -  2y) 


1  -  \/l  -  m2 

2|m! 


(2.14) 


& 

'A' 


the  maximum  value  for  P,  which  will  determine  the  stability  of  the  differencing  schemes,  is 


Pmax  — 


(1  -  \/l  -  /l2)Al 

2|HAy 


(2.15) 


Flows  with  larger  \jl  will  tend  to  be  less  stable. 

The  integration  was  performed  for  0  <  x  <  1  and  0  <  y  <  1.  Figure  2.1  shows  the  result 
of  the  integrations  when  Ax  =  Ay  =  0.05,  n  =  0.3  and  the  starting  values  were  given  by 
/(0,  y)  =  1  —  cos(?ry).  The  deviation  of  the  predicted  values  from  the  exact  values  are  also 
shown.  In  this  flow  the  velocity  component  Vv  is  both  small  and  varies  smoothly  with  y.  Both 
upwind  and  partial  upwind  differencing  are  stable  (Pma z  =  0.077)  but  the  instability  of  the 
central  differncing  scheme  is  not  clearly  manifested.  On  the  other  hand,  it  is  clear  that  the 
upwind  differencing  scheme  has  much  higher  truncation  error  relative  to  the  other  schemes. 

Figure  2.2  shows  the  results  of  integrations  with  Ax  =  0.1,  Ay  =  0.0333,  /r  =  0.9  and 
starting  values  given  by  the  pyramid  function  /( 0, y)  =  (1  —  [1  —  2y|)/2.  Again,  the  upwind  and 
partial  upwind  differencing  schemes  are  stable  (Pmaz  =  0.940),  but  in  this  case  the  instability 
of  the  central  differencing  scheme  is  clear.  The  rapid  change  in  the  y-derivative  of  Vv  causes 
oscillations  to  appear  in  the  central  difference  solution.  Once  more,  the  accuracy  of  the  partial 
upwind  scheme  is  the  best. 
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3  Stability  of  a  Hyperbolic  System  of 
Equations 


The  stability  of  a  system  of  first  order  equations  can  often  be  reduced  to  the  stability  of 
a  number  of  one  dimensional  equations.  Suppose  /„,n  =  1 ,N  are  variables  dependent  on 
x  and  y,  and  which  satisfy 

n=l  3 

where  A  and  B  are  N  x  N  matrices  and  A  is  non-singular.  Let  An  be  the  eigenvalues  and  v(n) 
be  N  linearly  independent  eigenvectors  of  A-1B:  i.e. 


£  £  'CM"1  =  *»•*’ 

;=1 t=l 


or  equivalently, 


E(A -  £*y)vjn)  =  0  (3.3) 

i=l 

If  all  N  eigenvalues  are  real  and  non-zero,  equation  (3.1)  is  said  to  be  hyperbolic.  If  they  are 
all  real  but  some  are  zero,  equation  (3.1)  is  said  to  be  mixed  parabolic-hyperbolic.  Henceforth 
each  An  will  be  assumed  to  be  real. 

The  dependent  variables,  /„,  can  be  expressed  as  a  linear  combination  of  the  eigenvectors 


Defining  vectors  w(n)  such  that 


one  can  represent  the  0k  explicitly: 


a  =  £  /-«£ 


Note,  too,  that 


*=i  t=i 


,»:ww 


as  can  be  verified  by  substituting  equation  (3.7)  into  equation  (3.2). 

After  mutiplication  by  A-1  and  substitution  for  /„,  equation  (3.1)  can  be  rewritten 

V"'  (n)d&n  (n)  .  dfin  _  /,  o') 

a7  +  v*  A"air  =  (3  8) 

n=l  9 

whence,  since  the  v(n)  are  linearly  independent, 

=  °  (3.9) 

Hence,  by  applying  a  stable  one-dimensional  differencing  scheme  to  each  of  the  /?„,  one  obtains 
a  stable  differencing  scheme  for  the  entire  linear  system.  The  differencing  scheme  for  /„  is 
recovered  using  equation  (3.4). 


N  S 


= 


(m),Jm)  ,  i±fl 


Y^D+fr  +  —^D~fr 


where 


D+fn  +  D-fn  „  D+fm-D-f„ 

9  L.  Unm-  n 


Dnm  =  £  vir)w^ar 


(3.10) 


(3.11) 


If  all  the  eigenvalues  are  of  the  same  sign,  it  is  possible  to  choose  all  the  am  to  be  equal 
and  have  a  stable  scheme.  This  has  the  advantage  that  the  cross-stream  differencing  is  the 
same  for  each  equation  in  the  system.  A  complete  upwind  scheme  would  have  am  =  sgn(Am) 
while  a  possible  extension  of  the  partial  upwind  scheme  of  the  previous  section  would  be  am  = 
max(Ai, . . . ,  A/v)Az/Ay.  By  using  the  maximum  eigenvalue  one  can  ensure  stability  by  choosing 
Ax  sufficiently  small  that  am  <  1.  When  all  the  am  are  equal,  equation  (3.10)  reduces  to  the 
simple  counterpart  of  equation  (2.3): 


(3.12) 


However,  if  the  eigenvalues  are  not  all  of  the  same  sign,  then  to  maintain  at  least  partial 
upwinding  to  provide  stability,  some  of  the  am  must  be  positive  and  some  must  be  negative. 
Thus,  in  general  it  is  not  possible  to  have  complete  stability  while  retaining  the  simplicity  of 
having  all  am  equal. 

The  logical  extension  of  the  partial  upwinding  scheme  of  the  previous  section,  is  to  set 
»m  =  Pm  Then 


(3.13) 
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whence 


to 
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D  =  A_1BAx/Ay  (3.14) 

The  truncation  error  when  D  is  defined  in  this  way  is  of  second  order  in  Ai  and  Ay. 


4  A  Differencing  Scheme  for  the 
Boundary  Layer  Equations 


As  mentioned  in  the  introduction,  the  boundary  layer  equations  discussed  here  are  assumed 
to  be  of  the  form  used  in  the  DREA  HLLFLO  programs[l].  These  are  a  system  of  three  coupled 
first  order  equations:  the  two  momentum  integral  equations  and  Head’s  entrainment  equation. 
They  can  be  written  in  the  form[2], 


where 


dQl 

dx 

dQ2 

-| — - —  =  production  terms 

da 

(4.1) 

QX1  dQX2 

— - 1 - - —  =  production  terms 

dx  da 

(4.2) 

dEl 

dx 

dE2  ,  . 

-i — - —  =  production  terms 

da 

(4.3) 

Q 1 

=  y/g6{ylhu  +  tv±lhu) 

(4.4) 

Q 2 

=  Vgf>{v2hn  +  tvx2hn) 

(4.5) 

qXI 

=  y/gStiV'hzi  +  tvxihzz) 

(4.6) 

QX2 

=  y/H?>t{V2hzi  +  tVX2hzz) 

(4.7) 

El 

=  y/gS{Vlkz  +  tvxikz) 

(4.8) 

E2 

=  yfis^kz+tv^kz) 

(4.9) 

6  is  the  boundary  layer  thickness,  and  t  is  the  tan  of  the  cross-flow  angle,  V1  and  V 2  are  the 
contravariant  components  of  the  potential  flow  velocity,  V11  and  V X2  are  the  contravariant 
components  of  a  vector  perpendicular  to  the  potential  flow  velocity  and  having  the  same  mag¬ 
nitude,  i  and  a  are  non-orthogonal  coordinates  which  run  over  the  surface  of  flow,  and  g  is  the 
determinant  of  the  metric  tensor  for  the  coordinates  ( x,s ).  The  values  of  V1,  V2,  VX1,  Vx2, 
and  g  are  known  for  any  given  z  and  a. 

The  profile  functions  hn,  h\z,  hzi,  hiz,  ki,  and  kz  depend  on  the  velocity  profile  assump¬ 
tions  made  for  the  boundary  layer  calculation.  In  general,  the  profile  functions  may  depend  on 
S,  t,  and  a  third  independent  variable  usually  taken  to  be  the  shape  factor,  JET,  though  more 
commonly  they  depend  only  upon  H  and  t  or  on  H  alone. 


A  simple  variant  of  upwind  differencing  can  be  used  to  integrate  the  boundary  layer 
equations  (4.1)  -  (4.3).  If  V2 /V1  and  (V2  +  tVX2)/[Vl  +  tVxl)  both  have  the  same  sign  (i.e. 
both  potential  streamlines  and  limiting  streamlines  lie  on  the  same  side  as  the  x  coordinate 
lines),  upwind  differencing  is  used  for  each  of  Q 2,  Q x2,  and  E2.  However,  if  they  have  different 
signs  (the  x  coordinate  line  bisects  the  two  streamlines),  no  unique  upwind  direction  is  identified: 
central  differencing  is  then  used.  Unfortunately,  as  has  been  demonstrated  above,  central 
differencing  is  unstable.  Moreover,  it  is  clear  from  Section  3  that  no  simple  scheme  in  which 
the  same  value  for  a  is  used  for  each  of  Q2,  Q x2,  and  E2  can  be  stable  when  the  x  coordinate 
line  lies  between  the  characteristic  directions  defined  by  Ai  and  A3. 

On  the  other  hand,  by  changing  independent  variables  to  the  boundary  layer  variables  S, 
H,  and  t,  the  boundary  layer  equations  may  be  rewritten  in  the  form  of  equation  (2.11)  with 
A  =  6,  A  =  H,  f3  =  t,  and 


(4.10) 


A  and  B  are  easy  to  determine  analytically  and,  since  A  is  only  3x3,  A-1  is  also  easily 
calculated.  Hence,  the  differencing  scheme  described  in  the  previous  section  may  be  calculated 
easily.  It  is  well-known  that  the  characteristics  of  the  boundary  layer  equations  lie  between  the 
potential  flow  streamlines  and  the  limiting  streamlines  at  the  hull  surface;  this  is  equivalent  to 
saying  that  the  eigenvalues  of  A-1B  lie  between  V2 /V1  and  (V2  +  tVr-L2)/(V1  +  tV"-11).  Hence, 
the  CFL  condition  is  satisfied  provided  that 
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dt 

dQ12 
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dt 
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dE2 

dS 

dH 

dt 

V2Ax  (V2  +  tVr-L2)Ax\  ^ 

V'Aa’  (V1  +  tVxl)As )  ~ 


(4.11) 


The  differencing  scheme  just  described  is  with  respect  to  the  independent  variables  6,  H 
and  t.  To  generate  corresponding  difference  schemes  for  Q2,  Qx2  and  E2,  the  correspondence 


dQ 2  dS_ 

da  da 

dQ 12  =B  d_H_ 
da  da 

dE2  dt_ 

da  .da  . 

is  used.  Equation  (3.10)  is  then  an  appropriate  difference  scheme  where 

ifi,f2,h)  =  {Q2,Q±2>E2) 

and 

D  =  B(A-1B)B-1Az/As  =  BA_1Az/As 


(4.12) 


(4.13) 


(4.14) 


If  it  is  assumed  that  the  profile  functions  depend  only  on  H ,  further  simplifications  are 
possible.  In  this  case  one  has 


where 


A  =  y/gT{V1C  +  tV'xlE)F  ;  B  =  y/gT(V2C  +  tF±2E)F 


10  0  10  0 

T  ='  0  t  0  ,  F  =  050 

0  0  1  0  0  5 


*12  **12 


(4.15) 


(4.16) 


c  —  /121  h'21  ti2i  >  E  —  hi2  h'22  2/122  (417) 

ki  k[  0  k2  k'2  ki 

The  primes  denote  differentiation  with  respect  to  H .  Substitution  into  equation  (4.15)  gives 


BA-1  =  TtV^C  +  tV'^EXV^C  +  tV^E)'-^-1 

=  T{V 2  +  tVi2EC-1)(V'1  +  tV'xlEC-l)-1T-1 


(4.18) 


Since  the  1  coordinate  lines  are  assumed  to  lie  roughly  along  the  potential  '"iow  streamlines, 
and  since  t  is  usually  (though  not  invariably)  significantly  smaller  than  1.0,  it  is  reasonable  to 
assume  that  tVxl  <<  V1  and  retain  only  first  order  terms  in  equation  (4.18). 


BA-1  =  T 


y 2  yly±2  _  y±ly2 

yi  **  (yi)2  T 


G  T' 


where 


G  =  EC-1 


(4.19) 


(4.20) 


and  U  is  the  magnitude  of  the  potential  flow  velocity.  The  last  equality  in  equation  (4.19) 
follows  from  the  relations  between  V  and  Vx  (see  Appendix  A. 5  in  Hally[2]). 

The  matrix  G  depends  only  on  H.  Its  components  may  be  written  explicitly  in  terms  of 
the  profile  functions  as  follows 


fci(fc^2i  ~  *2^21) 

-  hnk[) 
k2 

l+hTi 
_  huGn 

__  *1(^22^21  -  ^22^21)  -  ^22(^21^1  ~  h2ik\) 

J»2i(A'nfci  -  Anfcj) 


(4.21) 

(4.22) 

(4.23) 

(4.24) 
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G22 
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(4.25) 

h22(h'21hu 

-  h2ih'n)  -  hn{k'22h2i  -  h22h21) 

<J23 

h2i(h'nkl  -hnk'J 

(4.ZbJ 

G31 

=  Gn 

(4.27) 

G32 

=  G 12  -  1 

(4.28) 

G33 

=  G 13 

(4.29) 

where  use  has  been  made  of  the  identity 

h\2  =  h  21  +  k2  (4-30) 


For  the  particular  case  of  power  law  profiles  with  Mager’s  cross-flow  assumption  (see 
Hally  [2]),  the  profile  functions  are 
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and  the  components  of  G  are 
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(4.33) 

(4.34) 
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(4.36) 
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(4.38) 
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(4.43) 
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To  a  very  good  approximation,  G  is  linear  in  H  over  the  range  of  interest  ( H  is  usually 
in  the  range  1.2  to  1.5).  A  linear  approximation  for  G  was  generated  using  a  least  squares  fit 
for  the  range  H  €  [1.2, 1.5],  yielding 


(4.44) 


The  components  of  G  and  the  approximations  to  them  are  shown  in  Figure  4.1. 

Substitution  of  the  above  approximations  into  equation  (3.10)  with  the  identifications  in 
equation  (4.13)  yields  the  following  approximation  to  the  stable  scheme  with  an  —  XnAx/As; 
therefore  it  should  be  stable  for  all  reasonable  values  of  the  boundary  layer  parameters. 
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(4.45) 
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(4.46) 


(4.47) 


This  scheme  is  efficient  to  calculate  (since  all  the  boundary  layer  parameters  are  known), 
has  small  truncation  error  and  small  artifical  viscosity,  and  is  continuous  with  respect  to  all 
boundary  layer  and  potential  flow  parameters. 

When  the  cross-flow  is  large,  when  the  shape  factor  lies  outside  the  range  [1.2, 1.5],  or 
when  profiles  other  than  the  power  law  profiles  are  used,  a  stable  differencing  scheme  can  be 
generating  by  calculating  BA-1  explicitly  and  using  equations  (3.10)  and  (4.14).  While  not  as 
efficient  as  using  equations  (4.45)  -  (4.47),  the  computational  effort  is  not  high  as  A  and  B  are 
only  3x3. 


5  Concluding  Remarks 


A  differencing  scheme  for  the  cross-stream  derivatives  of  integral  boundary  layer  equations 
has  been  discussed.  In  comparison  with  upwinding  differencing,  the  new  scheme  has  better 
stability  properties,  smaller  truncation  error,  smaller  artificial  viscosity,  and  is  continuous  with 
respect  to  the  direction  of  flow.  Hence,  it  is  a  more  accurate  and  more  reliable  method  for  the 
solution  of  the  boundary  layer  equations. 

When  applied  to  the  boundary  layer  equations  used  by  the  DREA  HLLFLO  programs, 
it  has  been  shown  that  the  improved  differencing  scheme  can  be  written  in  a  form  which  is 
straightforward  and  efficient  to  calculate. 
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